Integration of multimodal cues does not alter mean but reduces variance in avian responses to predators: a systematic review and meta-analysis

1 Setting up

Code

# install.packages("pacman")
# pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)

# devtools::install_github("daniel1noble/orchaRd", force = TRUE)
library(pacman)

p_load(
    tidyverse,
    here,
    lme4,
    orchaRd,
    metafor,
    patchwork,
    alluvial,
    ggalluvial,
    ape,
    easyalluvial,
    clubSandwich,
    emmeans,
    MuMIn,
    kableExtra,
    pander
)

# library(tidyverse)
# library(here)
# library(lme4)
# library(orchaRd)
# library(metafor)
# library(patchwork)
# library(alluvial)
# library(ggalluvial)
# library(easyalluvial)
# library(ape)
# library(clubSandwich)
# library(emmeans)
# library(MuMIn)
# library(kableExtra)
# library(pander)
# making metafor talk to MuMIn
eval(metafor:::.MuMIn)
# install.packages("pak")
# pak::pak("MichelNivard/gptstudio")

2 Getting data loaded

Code

#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))
#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))
dat_full <- read.csv(here("data/dat_19_07_2023.csv"))

# add phylogenetic tree - only topologies
# TODO? - we could get better tree from birdtree.org
# we can do 50 different trees as in 
# https://academic.oup.com/sysbio/article/68/4/632/5267840
tree_top <- read.tree(here("R/birds_MA.tre"))

# tree with branch lengths
tree <- compute.brlen(tree_top)
#plot(tree)
# turning it into a correlation matrix
cor_tree <- vcv(tree,corr=T)

3 Custom functions

Code

# custom functions

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name",
                                       "\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "homo") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  VCV1 <- vcalc(vi = dat2$Vd,
             cluster = dat2$SubjectID,
             rho = 0.5)
  #model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)
  # meta-regression: contrasts 
  # helper function to run metafor meta-regression
  run_rma1 <- function(name) {
    rma.mv(yi = SMD, 
         V = VCV1, 
         mods = ~ relevel(dat2[[mod]], ref = name), 
         random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat2)
   }

    run_rma2 <- function(name) {
    rma.mv(yi = SMD,
         V = VCV1,
         mods = ~ relevel(dat2[[mod]], ref = name),
         random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             formula(paste("~", mod, "| Obs_ID"))),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat2)
   }

  # results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
  # this does not work for hetero model?
  if (type == "homo"){
  model_all <- map(level_names[-length(level_names)], run_rma1)
  } else {
  model_all <- map(level_names[-length(level_names)], run_rma2)
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# test 
#all_models(mod1, mod = "Treat_mod")
#all_models(mod1s, mod = "Treat_mod")
#all_models(mod1s2, mod = "Treat_mod")

4 Preparing data set

Code

# function for calculating variance
Vd_func <- function(d, n1, n2, design, r = 0.5){
  # independent design
  if(design == "among"){
    var <- (n1 + n2) / (n1*n2) + d^2 / (2 * (n1 + n2 - 2)) # variance
  } else { # dependent design
    var <- 2*(1-r) / n1 + d^2 / (2*(n1 - 1)) # variance
  }
  var # return variance
}

# getting Hedges' g - get small size bias corrected effect size
dat_full$SMD <- dat_full$d / (1 - 3/(4 * (dat_full$NTreat + dat_full$Ncontrol) - 9))

# flipping d 
dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection


# calucating Vd
dat_full$Vd <- with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))


# extra useful function
# function for getting mean and sd from median, quartiles and sample size
# get_mean_sd <- function(median, q1, q3, n){
#   sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd
#   mean <- (median + q1 + q3)/3 # mean
#   c(mean, sd)
# }


# observation id
dat_full$Obs_ID <- 1:nrow(dat_full)
dat_full$Phylo <- gsub(" ", "_", dat_full$FocalSpL)

# filtering very large variance and also very small sample size
dat_int <- dat_full %>% filter(Vd < 10 & Ncontrol > 2 & NTreat > 2)

#dim(dat_full)
#dim(dat_int)


# sorting out modality stuff
# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) 

dat_int %>% mutate(Treat_mod = case_when(Treatment == "A" ~ "A",
                                          Treatment == "AV" ~ "AV",
                                          Treatment == "AVG" ~ "AV",
                                          Treatment == "AVM" ~ "AV",
                                          Treatment == "L" ~ "AVO",
                                          Treatment == "O" ~ "O",
                                          Treatment == "OV" ~ "OV",
                                          Treatment == "V" ~ "V",
                                          Treatment == "VG" ~ "V",
                                          Treatment == "VM" ~ "V",
                                          Treatment == "VP" ~ "V"),
                    # into how many
                    Treat_No = case_when(Treatment == "A" ~ 1,
                                         Treatment == "AV" ~ 2,
                                         Treatment == "AVG" ~ 2,
                                         Treatment == "AVM" ~ 2,
                                         Treatment == "L" ~ 3,
                                         Treatment == "O" ~ 1,
                                         Treatment == "OV" ~ 2,
                                         Treatment == "V" ~ 1,
                                         Treatment == "VG" ~ 1,
                                         Treatment == "VM" ~ 1,
                                         Treatment == "VP" ~ 1),
                    # des it have some add-ons
                    Add_on = case_when(Treatment == "A" ~ "No",
                                         Treatment == "AV" ~ "No",
                                         Treatment == "AVG" ~ "Yes",
                                         Treatment == "AVM" ~ "Yes",
                                         Treatment == "L" ~ "No",
                                         Treatment == "O" ~ "No",
                                         Treatment == "OV" ~ "No",
                                         Treatment == "V" ~ "No",
                                         Treatment == "VG" ~ "Yes",
                                         Treatment == "VM" ~ "Yes",
                                         Treatment == "VP" ~ "Yes"),

                      ) -> dat

# creating data just for A, V, and AV 
dat_short <- dat %>% filter(Treat_mod == "A" | Treat_mod == "V" | Treat_mod == "AV")

# for add-on, we only need V and AV
dat_short_add <- dat %>% filter(Treat_mod == "AV" | Treat_mod == "V")


dat <- dat %>%
  mutate_if(is.character, as.factor)

# reordering factors for better visualization

# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))

5 Exploratory visualization

For Treat_mod (Treatment), we will only visualise A, V, and AV as O $r $, OV, and AVO are much rarer. But for Type (Trait type), we will use all data.

Treatment vs. Trait types

Code

# reordering dat_shrot$Treat_mod for better visualisation
dat_short$Treat_mod <- factor(dat_short$Treat_mod, levels = c("A", "V", "AV"))

# Treatment vs Type
dat_short %>% group_by(Treat_mod, Type) %>%
  summarise(n = n()) -> tab

#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Type)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and trait type")

Treatment vs. Treatment duration

Code

# tuning Treatment duration into a binary variable
dat_short %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat_short

# reformatting data
dat_short %>% group_by(Treat_mod, TDration) %>%
  summarise(n = n()) -> tab

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = TDration)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and duration of treatment")

Treatment vs. Sex

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Sex) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Predator type

Code

# Treat_mod vs Design
dat_short %>% filter(PredTo != "") %>%
  group_by(Treat_mod, PredTo) %>%
  summarise(n = n()) -> tab

tab$PredTo <- factor(tab$PredTo, levels = c("A", "N", "B"), labels = c("Adult", "Nestling", "Both"))

#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = PredTo)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Design

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Design) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Design)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Season

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Season) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Season)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Setting

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Setting) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Setting)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs Control type

Code

# Treat_mod vs ControlType

dat_short %>% filter(ControlType != "mix") %>% group_by(Treat_mod, ControlType) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = ControlType)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 6, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and control type")

Trait type vs. Treatment duration

Code

# Type vs duration_days

# turn duration_days into a binary factor: less than 1 and more than 1

dat %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat


dat %>% group_by(Type, TDration) %>%
  summarise(n = n()) -> tab2

# using ggaruvial
ggplot(tab2,
       aes(y = n,
           axis1 = Type,
           axis2 = TDration)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and duration of treatment")

Trait type vs. Sex

Code

# Type vs duration_days

# turn duration_days into a binary factor: less than 1 and more than 1

dat %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat

dat %>% group_by(Type, Sex) %>%
  summarise(n = n()) -> tab

# reordering
tab$Sex <- factor(tab$Sex, levels = c("F", "M", "both"), labels = c("Female", "Male", "Both"))


# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Type,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and duration of treatment")

6 Meta-analysis

All random effects (one tree)

Code

##| warning: false
# VCV matrix
# TODO - new phylogeny

VCV <- vcalc(vi = dat$Vd,
             cluster = dat$SubjectID,
             rho = 0.5)

mod_ma <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(#~1 | Phylo,
                     ~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | SubjectID, # incoprated as VCV
                     ~1 | Obs_ID),
      # R = list(Phylo = cor_tree),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod_ma)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-894.6453  1789.2906  1799.2906  1821.5901  1799.3853   

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0111  0.1053     90     no   FocalSpL 
sigma^2.2  0.1450  0.3808    116     no      RecNo 
sigma^2.3  0.0000  0.0001    153     no  SubjectID 
sigma^2.4  0.6872  0.8290    640     no     Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4139  0.0653  6.3401  639  <.0001  0.2857  0.5421  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# TODO - think about whether we add this or not
robust(mod_ma, cluster = dat$SubjectID)

Multivariate Meta-Analysis Model (k = 640; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0111  0.1053     90     no   FocalSpL 
sigma^2.2  0.1450  0.3808    116     no      RecNo 
sigma^2.3  0.0000  0.0001    153     no  SubjectID 
sigma^2.4  0.6872  0.8290    640     no     Obs_ID 

Number of estimates:   640
Number of clusters:    153
Estimates per cluster: 1-29 (mean: 4.18, median: 2)

Model Results:

estimate      se¹    tval¹   df¹    pval¹   ci.lb¹   ci.ub¹      
  0.4139  0.0540   7.6691   152   <.0001   0.3073   0.5206   *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1) results based on cluster-robust inference (var-cov estimator: CR1,
   approx t-test and confidence interval, df: residual method)

Code

round(i2_ml(mod_ma), 2)

    I2_Total  I2_FocalSpL     I2_RecNo I2_SubjectID    I2_Obs_ID 
       92.81         1.22        15.96         0.00        75.63 

Code

orchard_plot(mod_ma,
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

All random effects (100 trees)

Code

# TODO

Reduced model

Code

# reduced model

mod_ma2 <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | Obs_ID),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod_ma2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-894.6453  1789.2906  1797.2906  1815.1302  1797.3537   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0111  0.1052     90     no  FocalSpL 
sigma^2.2  0.1450  0.3808    116     no     RecNo 
sigma^2.3  0.6872  0.8290    640     no    Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4139  0.0653  6.3401  639  <.0001  0.2857  0.5421  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(i2_ml(mod_ma2), 2)

   I2_Total I2_FocalSpL    I2_RecNo   I2_Obs_ID 
      92.81        1.22       15.96       75.63 

Code

# comparing two models
#anova(mod0, mod0r)

7 Meta-regression: uni-moderator

Treatmeant with all data

Code

mod_trt <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Treat_mod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod_trt)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-887.7086  1775.4172  1793.4172  1833.4856  1793.7056   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0074  0.0862     90     no  FocalSpL 
sigma^2.2  0.1563  0.3954    116     no     RecNo 
sigma^2.3  0.6907  0.8311    640     no    Obs_ID 

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 634) = 6.8763, p-val < .0001

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub      
Treat_modA      0.4019  0.0971  4.1381  634  <.0001   0.2112  0.5926  *** 
Treat_modAV     0.4723  0.1345  3.5128  634  0.0005   0.2083  0.7363  *** 
Treat_modAVO    0.5955  0.4487  1.3270  634  0.1850  -0.2857  1.4767      
Treat_modO      0.2577  0.2930  0.8797  634  0.3793  -0.3176  0.8330      
Treat_modOV     0.0185  0.3829  0.0482  634  0.9615  -0.7334  0.7703      
Treat_modV      0.4324  0.1043  4.1442  634  <.0001   0.2275  0.6373  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_trt)*100, 2)

   R2_marginal R2_conditional 
          0.64          19.68 

Code

orchard_plot(mod_trt, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table
all_models(mod_trt, mod = "Treat_mod")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.402	0.211	0.593	0.000	-1.423	2.227
AV	0.472	0.208	0.736	0.000	-1.362	2.307
AVO	0.595	-0.286	1.477	0.185	-1.422	2.613
O	0.258	-0.318	0.833	0.379	-1.646	2.162
OV	0.018	-0.733	0.770	0.962	-1.946	1.983
V	0.432	0.228	0.637	0.000	-1.394	2.259
A-AV	0.070	-0.248	0.389	0.664	-1.772	1.913
A-AVO	0.194	-0.706	1.093	0.673	-1.832	2.219
A-O	-0.144	-0.748	0.460	0.639	-2.057	1.769
A-OV	-0.383	-1.158	0.392	0.332	-2.357	1.590
A-V	0.031	-0.240	0.301	0.825	-1.805	1.866
AV-AVO	0.123	-0.797	1.043	0.793	-1.912	2.158
AV-O	-0.215	-0.847	0.418	0.506	-2.137	1.708
AV-OV	-0.454	-1.250	0.343	0.264	-2.436	1.528
AV-V	-0.040	-0.340	0.260	0.794	-1.880	1.800
AVO-O	-0.338	-1.390	0.714	0.529	-2.436	1.760
AVO-OV	-0.577	-1.711	0.557	0.318	-2.717	1.563
AVO-V	-0.163	-1.067	0.741	0.723	-2.191	1.865
O-OV	-0.239	-1.186	0.707	0.620	-2.286	1.808
O-V	0.175	-0.434	0.784	0.574	-1.740	2.089
OV-V	0.414	-0.364	1.192	0.296	-1.561	2.389

Treatment with Audio and Visual data

Code

#homoscedastic model

VCVs <- vcalc(vi = dat_short$Vd,
             cluster = dat_short$SubjectID,
             rho = 0.5)


mod_trt_1 <- rma.mv(yi = SMD,
                   V = VCVs,
                   random = list(~1|FocalSpL,
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mod = ~ Treat_mod - 1, 
                   test = "t",
                   method = "REML", 
                   sparse = TRUE,
                   data = dat_short)

summary(mod_trt_1)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-847.7054  1695.4108  1707.4108  1733.7623  1707.5532   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0087  0.0932     81     no  FocalSpL 
sigma^2.2  0.1700  0.4123    104     no     RecNo 
sigma^2.3  0.7138  0.8449    600     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 597) = 12.6006, p-val < .0001

Model Results:

             estimate      se    tval   df    pval   ci.lb   ci.ub      
Treat_modA     0.4055  0.0996  4.0715  597  <.0001  0.2099  0.6012  *** 
Treat_modV     0.4386  0.1069  4.1041  597  <.0001  0.2287  0.6485  *** 
Treat_modAV    0.4782  0.1379  3.4687  597  0.0006  0.2074  0.7490  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_trt_1)*100, 2)

   R2_marginal R2_conditional 
          0.08          20.09 

Code

orchard_plot(mod_trt_1, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table (get all the constrasts)
all_models(mod_trt_1, mod = "Treat_mod")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.406	0.210	0.601	0.000	-1.460	2.271
V	0.439	0.229	0.648	0.000	-1.429	2.306
AV	0.478	0.207	0.749	0.001	-1.397	2.353
A-V	0.033	-0.244	0.310	0.815	-1.843	1.909
A-AV	0.073	-0.254	0.399	0.662	-1.811	1.957
V-AV	0.040	-0.267	0.346	0.800	-1.841	1.920

Code

# modeling heteroscedasticity
mod_trt_2 <- rma.mv(yi = SMD, 
                V = VCVs, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ Treat_mod | Obs_ID), 
                mod = ~ Treat_mod - 1, 
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short)

summary(mod_trt_2)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-838.2594  1676.5188  1692.5188  1727.6541  1692.7637   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0086  0.0927     81     no  FocalSpL 
sigma^2.2  0.1764  0.4200    104     no     RecNo 

outer factor: Obs_ID    (nlvls = 600)
inner factor: Treat_mod (nlvls = 3)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.7479  0.8648    302     no      A 
tau^2.2    0.8583  0.9264    190     no      V 
tau^2.3    0.3460  0.5882    108     no     AV 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 597) = 12.4198, p-val < .0001

Model Results:

             estimate      se    tval   df    pval   ci.lb   ci.ub      
Treat_modA     0.4100  0.1013  4.0495  597  <.0001  0.2112  0.6089  *** 
Treat_modV     0.4400  0.1108  3.9729  597  <.0001  0.2225  0.6576  *** 
Treat_modAV    0.4482  0.1229  3.6459  597  0.0003  0.2068  0.6897  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod_trt_2, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table (get all the constrasts)
all_models(mod_trt_2, mod = "Treat_mod")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.410	0.211	0.609	0.000	-1.497	2.317
V	0.440	0.223	0.658	0.000	-1.578	2.458
AV	0.448	0.207	0.690	0.000	-1.003	1.900
A-V	0.033	-0.244	0.310	0.815	-1.843	1.909
A-AV	0.073	-0.254	0.399	0.662	-1.811	1.957
V-AV	0.040	-0.267	0.346	0.800	-1.841	1.920

Code

# comparison models
anova(mod_trt_1, mod_trt_2) # modeling heteroscedasticity is better 

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1692.5188 1727.6541 1692.7637 -838.2594                NA 
Reduced  6 1707.4108 1733.7623 1707.5532 -847.7054 18.8921 <.0001 NA 

Treatment with additions

Code

# the effect of additions
# this is a part of sensitivity analysis

dat_short_add$Treat_add <- as.factor(paste(dat_short_add$Treat_mod , 
                                 dat_short_add$Add_on, 
                                 sep = "-"))

VCVs2 <- vcalc(vi = dat_short_add$Vd,
             cluster = dat_short_add$SubjectID,
             rho = 0.5)

mod_add <- rma.mv(yi = SMD, 
                V = VCVs2, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ 1 | Obs_ID), 
                mod = ~ Treat_add - 1,
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short_add)

summary(mod_add) 

Multivariate Meta-Analysis Model (k = 298; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-422.9579   845.9158   859.9158   885.7008   860.3074   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0000     39     no  FocalSpL 
sigma^2.2  0.1432  0.3784     58     no     RecNo 
sigma^2.3  0.6666  0.8165    298     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 294) = 6.2038, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
Treat_addAV-No     0.4270  0.1410  3.0289  294  0.0027   0.1495  0.7045   ** 
Treat_addAV-Yes    0.5716  0.2828  2.0210  294  0.0442   0.0150  1.1283    * 
Treat_addV-No      0.4262  0.1095  3.8912  294  0.0001   0.2106  0.6418  *** 
Treat_addV-Yes     0.2930  0.2091  1.4017  294  0.1621  -0.1184  0.7045      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod_add, 
             mod = "Treat_add",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table (get all the constrasts)
all_models(mod_add, mod = "Treat_add")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
AV-No	0.427	0.150	0.704	0.003	-1.366	2.220
AV-Yes	0.572	0.015	1.128	0.044	-1.285	2.428
V-No	0.426	0.211	0.642	0.000	-1.358	2.210
V-Yes	0.293	-0.118	0.704	0.162	-1.525	2.111
AV-No-AV-Yes	0.145	-0.466	0.755	0.642	-1.729	2.018
AV-No-V-No	-0.001	-0.324	0.322	0.996	-1.801	1.800
AV-No-V-Yes	-0.134	-0.619	0.351	0.587	-1.970	1.702
AV-Yes-V-No	-0.145	-0.734	0.443	0.627	-2.012	1.721
AV-Yes-V-Yes	-0.279	-0.936	0.379	0.405	-2.168	1.611
V-No-V-Yes	-0.133	-0.578	0.312	0.556	-1.959	1.693

Treatment as an ordinal variable

Code

# testing the number of stimuli

mod_ord <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Treat_No, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod_ord)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.4118  1786.8236  1796.8236  1819.1153  1796.9185   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0117  0.1081     90     no  FocalSpL 
sigma^2.2  0.1469  0.3833    116     no     RecNo 
sigma^2.3  0.6875  0.8292    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 0.0952, p-val = 0.7577

Model Results:

          estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt     0.3696  0.1603  2.3057  638  0.0214   0.0548  0.6843  * 
Treat_No    0.0366  0.1185  0.3086  638  0.7577  -0.1962  0.2694    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod_ord,
             mod = "Treat_No",
             group = "RecNo",
             xlab = "The number of simuli",
             g = TRUE)

Trait type

Code

#homoscedastic model
# Type of responses 
mod_type1 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Type - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod_type1)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-886.7295  1773.4590  1785.4590  1812.1996  1785.5923   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     90     no  FocalSpL 
sigma^2.2  0.1322  0.3635    116     no     RecNo 
sigma^2.3  0.6902  0.8308    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 18.0966, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.4866  0.0677  7.1917  637  <.0001   0.3537  0.6194  *** 
TypeLifeHistory    0.3089  0.1176  2.6263  637  0.0088   0.0779  0.5398   ** 
TypePhysiology     0.0284  0.1285  0.2207  637  0.8254  -0.2240  0.2807      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_type1)*100, 2)

   R2_marginal R2_conditional 
          3.15          18.71 

Code

orchard_plot(mod_type1,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table
all_models(mod_type1, mod = "Type")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Behaviour	0.487	0.354	0.619	0.000	-1.299	2.272
LifeHistory	0.309	0.078	0.540	0.009	-1.487	2.105
Physiology	0.028	-0.224	0.281	0.825	-1.770	1.827
Behaviour-LifeHistory	-0.178	-0.418	0.062	0.146	-1.975	1.619
Behaviour-Physiology	-0.458	-0.723	-0.193	0.001	-2.259	1.342
LifeHistory-Physiology	-0.280	-0.586	0.025	0.072	-2.087	1.526

Code

# heteroscadasticity model
mod_type2 <- rma.mv(yi = SMD, 
               V = VCV, 
               mod = ~ Type - 1, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod_type2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-870.7865  1741.5731  1757.5731  1793.2272  1757.8024   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0274  0.1656     90     no  FocalSpL 
sigma^2.2  0.1296  0.3600    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6577  0.8110    434     no    Behaviour 
tau^2.2    0.3176  0.5636    112     no  LifeHistory 
tau^2.3    1.2391  1.1131     94     no   Physiology 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 17.4552, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.5026  0.0720  6.9850  637  <.0001   0.3613  0.6439  *** 
TypeLifeHistory    0.3531  0.1047  3.3712  637  0.0008   0.1474  0.5587  *** 
TypePhysiology     0.0120  0.1549  0.0773  637  0.9384  -0.2922  0.3161      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod_type2, 
             mod = "Type",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table
all_models(mod_type2, mod = "Type")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Behaviour	0.503	0.361	0.644	0.000	-1.275	2.281
LifeHistory	0.353	0.147	0.559	0.001	-1.015	1.721
Physiology	0.012	-0.292	0.316	0.938	-2.328	2.352
Behaviour-LifeHistory	-0.178	-0.418	0.062	0.146	-1.975	1.619
Behaviour-Physiology	-0.458	-0.723	-0.193	0.001	-2.259	1.342
LifeHistory-Physiology	-0.280	-0.586	0.025	0.072	-2.087	1.526

Code

# heteroscadasticity model better than the homoscedasticity model
# note LifeHistory has lowest variation but this may be expected? 
# as it is less flexiable (e.g. the number of eggs?)
anova(mod_type1, mod_type2)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1757.5731 1793.2272 1757.8024 -870.7865                NA 
Reduced  6 1785.4590 1812.1996 1785.5923 -886.7295 31.8859 <.0001 NA 

Treatment duration

Code

# treatment duration

dat$ln_duration <- log(dat$duration_days)

mod_dur <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_dur)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-889.1410  1778.2821  1788.2821  1810.5738  1788.3770   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1518  0.3897    116     no     RecNo 
sigma^2.3  0.6888  0.8300    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 8.4146, p-val = 0.0039

Model Results:

             estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt        0.3710  0.0634   5.8526  638  <.0001   0.2465   0.4955  *** 
ln_duration   -0.0455  0.0157  -2.9008  638  0.0039  -0.0763  -0.0147   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_dur)*100, 2)

   R2_marginal R2_conditional 
          3.43          20.87 

Code

bubble_plot(mod_dur,
             mod = "ln_duration",
             group = "RecNo",
             xlab = "log(duration in days)",
             g = TRUE) +
    geom_point(data = dat,
    aes(x = ln_duration, y = SMD,
    color = Type,
    fill = Type,
    size = 1/sqrt(Vd)), alpha = 0.6) +
    scale_color_discrete() + 
    guides(color = "legend")

Code

#p + geom_point(aes(colour = Type))
#scale_colour_manual(values = c("red", "blue", "green"))

# extra  - looking at the interaction
mod_dur_int <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration*Type,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

# no significant interaction
summary(mod_dur_int)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.1624  1760.3248  1778.3248  1818.3933  1778.6133   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1427  0.3777    116     no     RecNo 
sigma^2.3  0.6882  0.8296    640     no    Obs_ID 

Test of Moderators (coefficients 2:6):
F(df1 = 5, df2 = 634) = 3.5544, p-val = 0.0035

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                        0.4394  0.0717   6.1286  634  <.0001   0.2986 
ln_duration                   -0.0411  0.0173  -2.3802  634  0.0176  -0.0750 
TypeLifeHistory               -0.3341  0.3591  -0.9302  634  0.3526  -1.0393 
TypePhysiology                -0.4638  0.1637  -2.8326  634  0.0048  -0.7853 
ln_duration:TypeLifeHistory    0.0837  0.1012   0.8265  634  0.4088  -0.1151 
ln_duration:TypePhysiology     0.0438  0.0495   0.8849  634  0.3765  -0.0534 
                               ci.ub      
intrcpt                       0.5802  *** 
ln_duration                  -0.0072    * 
TypeLifeHistory               0.3712      
TypePhysiology               -0.1423   ** 
ln_duration:TypeLifeHistory   0.2825      
ln_duration:TypePhysiology    0.1409      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

res<-mod_results(mod_dur_int, mod = "ln_duration", group = "RecNo", by = "Type")

bubble_plot(res,
             mod = "ln_duration",
             group = "RecNo",
             xlab = "log(duration in days)",
             condition.nrow = 3,
             g = TRUE) 

Sex

Code

# sex

mod_sex1 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Sex - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_sex1)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-890.5710  1781.1420  1793.1420  1819.8826  1793.2753   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0099  0.0993     90     no  FocalSpL 
sigma^2.2  0.1354  0.3680    116     no     RecNo 
sigma^2.3  0.6906  0.8310    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 15.0441, p-val < .0001

Model Results:

         estimate      se    tval   df    pval   ci.lb   ci.ub      
Sexboth    0.4706  0.0758  6.2110  637  <.0001  0.3218  0.6193  *** 
SexF       0.2579  0.1011  2.5501  637  0.0110  0.0593  0.4565    * 
SexM       0.4309  0.1325  3.2515  637  0.0012  0.1707  0.6912   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_sex1)*100, 2)

   R2_marginal R2_conditional 
          0.95          18.17 

Code

orchard_plot(mod_sex1,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_sex1, mod = "Sex")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Both	0.471	0.322	0.619	0.000	-1.331	2.272
F	0.258	0.059	0.456	0.011	-1.548	2.064
M	0.431	0.171	0.691	0.001	-1.383	2.245
Both-F	-0.213	-0.434	0.009	0.059	-2.022	1.596
Both-M	-0.040	-0.323	0.244	0.784	-1.857	1.778
F-M	0.173	-0.098	0.444	0.211	-1.643	1.989

Code

# heteroscedasticity

mod_sex2 <- rma.mv(yi = SMD, 
                V = VCV, 
                mod = ~ Sex - 1, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~Sex | Obs_ID), 
                struct = "DIAG",
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod_sex2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-867.9399  1735.8799  1751.8799  1787.5341  1752.1092   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0007     90     no  FocalSpL 
sigma^2.2  0.1454  0.3813    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Sex    (nlvls = 3)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.8865  0.9415    398     no   both 
tau^2.2    0.2883  0.5369    158     no      F 
tau^2.3    0.6740  0.8210     84     no      M 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 14.8707, p-val < .0001

Model Results:

         estimate      se    tval   df    pval   ci.lb   ci.ub      
Sexboth    0.4729  0.0781  6.0557  637  <.0001  0.3196  0.6263  *** 
SexF       0.2429  0.0834  2.9127  637  0.0037  0.0791  0.4066   ** 
SexM       0.4268  0.1278  3.3403  637  0.0009  0.1759  0.6777  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod_sex2,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_sex2, mod = "Sex", type = "hetero")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Both	0.473	0.320	0.626	0.000	-1.528	2.474
F	0.243	0.079	0.407	0.004	-1.061	1.546
M	0.427	0.176	0.678	0.001	-1.368	2.222
Both-F	-0.213	-0.434	0.009	0.059	-2.022	1.596
Both-M	-0.040	-0.323	0.244	0.784	-1.857	1.778
F-M	0.173	-0.098	0.444	0.211	-1.643	1.989

Code

anova(mod_sex1, mod_sex2)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1751.8799 1787.5341 1752.1092 -867.9399                NA 
Reduced  6 1793.1420 1819.8826 1793.2753 -890.5710 45.2621 <.0001 NA 

Trait categories

Code

# Category of responses


mod_cat <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Category - 1, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod_cat)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-879.5908  1759.1816  1781.1816  1830.1194  1781.6074   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1389  0.3727    116     no     RecNo 
sigma^2.3  0.6912  0.8314    640     no    Obs_ID 

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 632) = 7.2433, p-val < .0001

Model Results:

                          estimate      se    tval   df    pval    ci.lb 
CategoryAntiPredator        0.4772  0.0924  5.1669  632  <.0001   0.2959 
CategoryCondition           0.0027  0.1394  0.0191  632  0.9848  -0.2710 
CategoryCostlyBehaviours    0.4567  0.1784  2.5596  632  0.0107   0.1063 
CategoryHormones            0.1717  0.2914  0.5891  632  0.5560  -0.4006 
CategoryIntake              0.7554  0.1623  4.6539  632  <.0001   0.4366 
CategoryParentalCare        0.4103  0.1076  3.8133  632  0.0002   0.1990 
CategoryPhenology           0.1927  0.1799  1.0711  632  0.2845  -0.1606 
Categoryreproduction        0.3366  0.1287  2.6159  632  0.0091   0.0839 
                           ci.ub      
CategoryAntiPredator      0.6586  *** 
CategoryCondition         0.2763      
CategoryCostlyBehaviours  0.8071    * 
CategoryHormones          0.7439      
CategoryIntake            1.0741  *** 
CategoryParentalCare      0.6216  *** 
CategoryPhenology         0.5459      
Categoryreproduction      0.5892   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_cat)*100, 2)

   R2_marginal R2_conditional 
          4.12          20.17 

Code

orchard_plot(mod_cat, 
             mod = "Category",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             angle = 45,
             branch.size = 3)

Code

# result table
all_models(mod_cat, mod = "Category")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
AntiPredator	0.477	0.296	0.659	0.000	-1.321	2.276
Condition	0.003	-0.271	0.276	0.985	-1.807	1.813
CostlyBehaviours	0.457	0.106	0.807	0.011	-1.366	2.280
Hormones	0.172	-0.401	0.744	0.556	-1.707	2.050
Intake	0.755	0.437	1.074	0.000	-1.062	2.573
ParentalCare	0.410	0.199	0.622	0.000	-1.391	2.212
Phenology	0.193	-0.161	0.546	0.285	-1.631	2.016
Reproduction	0.337	0.084	0.589	0.009	-1.470	2.143
AntiPredator-Condition	-0.475	-0.788	-0.161	0.003	-2.291	1.342
AntiPredator-CostlyBehaviours	-0.021	-0.389	0.348	0.913	-1.847	1.806
AntiPredator-Hormones	-0.306	-0.901	0.290	0.314	-2.191	1.580
AntiPredator-Intake	0.278	-0.063	0.619	0.110	-1.543	2.100
AntiPredator-ParentalCare	-0.067	-0.329	0.195	0.616	-1.875	1.741
AntiPredator-Phenology	-0.285	-0.662	0.093	0.140	-2.113	1.544
AntiPredator-Reproduction	-0.141	-0.428	0.146	0.336	-1.953	1.671
Condition-CostlyBehaviours	0.454	0.018	0.891	0.042	-1.388	2.296
Condition-Hormones	0.169	-0.445	0.783	0.589	-1.723	2.061
Condition-Intake	0.753	0.353	1.152	0.000	-1.081	2.586
Condition-ParentalCare	0.408	0.086	0.730	0.013	-1.410	2.226
Condition-Phenology	0.190	-0.228	0.608	0.373	-1.647	2.027
Condition-Reproduction	0.334	-0.004	0.671	0.053	-1.487	2.155
CostlyBehaviours-Hormones	-0.285	-0.952	0.382	0.402	-2.195	1.625
CostlyBehaviours-Intake	0.299	-0.160	0.757	0.202	-1.548	2.146
CostlyBehaviours-ParentalCare	-0.046	-0.448	0.356	0.821	-1.880	1.787
CostlyBehaviours-Phenology	-0.264	-0.753	0.225	0.289	-2.119	1.591
CostlyBehaviours-Reproduction	-0.120	-0.543	0.303	0.578	-1.959	1.718
Hormones-Intake	0.584	-0.068	1.235	0.079	-1.320	2.488
Hormones-ParentalCare	0.239	-0.362	0.839	0.435	-1.649	2.126
Hormones-Phenology	0.021	-0.625	0.667	0.949	-1.881	1.923
Hormones-Reproduction	0.165	-0.447	0.777	0.597	-1.726	2.056
Intake-ParentalCare	-0.345	-0.720	0.030	0.071	-2.173	1.483
Intake-Phenology	-0.563	-1.031	-0.094	0.019	-2.412	1.287
Intake-Reproduction	-0.419	-0.817	-0.020	0.039	-2.252	1.414
ParentalCare-Phenology	-0.218	-0.611	0.176	0.278	-2.050	1.614
ParentalCare-Reproduction	-0.074	-0.374	0.227	0.630	-1.888	1.740
Phenology-Reproduction	0.144	-0.220	0.508	0.438	-1.682	1.970

Predactor guild

Code

# Redactor guild
# Fish only N = 1 so remove and also remove NS (not specified)

dat %>% filter(PredGuild != "Fish" & PredGuild != "NS") ->
  data_preg

VCVs3 <- vcalc(vi = data_preg$Vd,
             cluster = data_preg$SubjectID,
             rho = 0.5)

mod_preg <- rma.mv(yi = SMD, 
               V = VCVs3, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredGuild - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = data_preg)

summary(mod_preg)

Multivariate Meta-Analysis Model (k = 620; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-857.4707  1714.9413  1728.9413  1759.9040  1729.1255   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0007  0.0271     85     no  FocalSpL 
sigma^2.2  0.1517  0.3895    110     no     RecNo 
sigma^2.3  0.6839  0.8270    620     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 616) = 10.0849, p-val < .0001

Model Results:

                  estimate      se    tval   df    pval    ci.lb   ci.ub      
PredGuildBird       0.4520  0.0732  6.1774  616  <.0001   0.3083  0.5957  *** 
PredGuildMammal     0.2681  0.1279  2.0957  616  0.0365   0.0169  0.5193    * 
PredGuildmixed      0.1120  0.2485  0.4508  616  0.6523  -0.3760  0.6001      
PredGuildReptile    0.3220  0.2800  1.1501  616  0.2506  -0.2279  0.8719      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_preg)*100, 2)

   R2_marginal R2_conditional 
          1.70          19.62 

Code

# R2_marginal R2_conditional 
#           1.70          19.62 

orchard_plot(mod_preg, 
             mod = "PredGuild",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_preg, mod = "PredGuild")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Bird	0.452	0.308	0.596	0.000	-1.350	2.254
Mammal	0.268	0.017	0.519	0.037	-1.545	2.082
Mixed	0.112	-0.376	0.600	0.652	-1.749	1.973
Reptile	0.322	-0.228	0.872	0.251	-1.556	2.200
Bird-Mammal	-0.184	-0.451	0.083	0.176	-2.000	1.632
Bird-Mixed	-0.340	-0.847	0.167	0.189	-2.206	1.526
Bird-Reptile	-0.130	-0.695	0.435	0.651	-2.013	1.753
Mammal-Mixed	-0.156	-0.705	0.392	0.577	-2.034	1.722
Mammal-Reptile	0.054	-0.550	0.657	0.861	-1.841	1.949
Mixed-Reptile	0.210	-0.525	0.945	0.575	-1.731	2.151

Setting

Code

# Setting

mod_set <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Setting - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod_set)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-891.0585  1782.1169  1794.1169  1820.8575  1794.2502   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0172  0.1311     90     no  FocalSpL 
sigma^2.2  0.1378  0.3712    116     no     RecNo 
sigma^2.3  0.6868  0.8288    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 14.2658, p-val < .0001

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub 
SettingField           0.4245  0.0721   5.8846  637  <.0001   0.2829  0.5662 
SettingLab             0.4764  0.1575   3.0247  637  0.0026   0.1671  0.7857 
SettingSemi-natural   -0.2054  0.4204  -0.4886  637  0.6253  -1.0308  0.6201 
                         
SettingField         *** 
SettingLab            ** 
SettingSemi-natural      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_set)*100, 2)

   R2_marginal R2_conditional 
          0.70          18.98 

Code

# R2_marginal R2_conditional 
#           0.70          18.98 

orchard_plot(mod_set, 
             mod = "Setting",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_set, mod = "Setting")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Field	0.425	0.283	0.566	0.000	-1.383	2.232
Lab	0.476	0.167	0.786	0.003	-1.352	2.304
Semi-natural	-0.205	-1.031	0.620	0.625	-2.187	1.776
Field-Lab	0.052	-0.283	0.387	0.761	-1.781	1.884
Field-Semi-natural	-0.630	-1.463	0.204	0.138	-2.615	1.355
Lab-Semi-natural	-0.682	-1.559	0.195	0.127	-2.686	1.322

Season

Code

# Season
mod_4 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Season - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)


summary(mod_4)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.1745  1784.3489  1794.3489  1816.6406  1794.4439   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0187  0.1368     90     no  FocalSpL 
sigma^2.2  0.1390  0.3728    116     no     RecNo 
sigma^2.3  0.6875  0.8291    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 21.2697, p-val < .0001

Model Results:

                    estimate      se    tval   df    pval   ci.lb   ci.ub      
Seasonbreeding        0.3589  0.0783  4.5831  638  <.0001  0.2051  0.5127  *** 
Seasonnon-breeding    0.5596  0.1123  4.9845  638  <.0001  0.3392  0.7801  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_4)*100, 2)

   R2_marginal R2_conditional 
          0.95          19.44 

Code

orchard_plot(mod_4,
             mod = "Season",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_4, mod = "Season")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Breeding	0.359	0.205	0.513	0.000	-1.453	2.171
Non-breeding	0.560	0.339	0.780	0.000	-1.259	2.378
Breeding-Non-breeding	0.201	-0.058	0.460	0.128	-1.623	2.025

Design

Code

# Design
mod_d <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Design - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_d)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.3891  1786.7783  1796.7783  1819.0700  1796.8732   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1481  0.3848    116     no     RecNo 
sigma^2.3  0.6879  0.8294    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 19.9969, p-val < .0001

Model Results:

              estimate      se    tval   df    pval   ci.lb   ci.ub      
Designamong     0.4061  0.0874  4.6483  638  <.0001  0.2345  0.5777  *** 
Designwithin    0.4197  0.0839  5.0006  638  <.0001  0.2549  0.5844  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_d)*100, 2)

   R2_marginal R2_conditional 
          0.01          18.72 

Code

 # R2_marginal R2_conditional 
 #          0.01          18.72 

orchard_plot(mod_d,
             mod = "Design",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_d, mod = "Design")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Among	0.406	0.235	0.578	0.000	-1.409	2.221
Within	0.420	0.255	0.584	0.000	-1.394	2.234
Among-Within	0.014	-0.204	0.231	0.903	-1.806	1.833

Response period

Code

# Response period
mod_resp <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ResponsePeriod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_resp)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.5749  1785.1498  1795.1498  1817.4415  1795.2448   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0138  0.1176     90     no  FocalSpL 
sigma^2.2  0.1435  0.3788    116     no     RecNo 
sigma^2.3  0.6883  0.8296    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 20.8414, p-val < .0001

Model Results:

                      estimate      se    tval   df    pval   ci.lb   ci.ub 
ResponsePeriodafter     0.3500  0.0851  4.1108  638  <.0001  0.1828  0.5171 
ResponsePeriodduring    0.4938  0.0903  5.4708  638  <.0001  0.3166  0.6711 
                          
ResponsePeriodafter   *** 
ResponsePeriodduring  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_resp)*100, 2)

   R2_marginal R2_conditional 
          0.61          19.10 

Code

#   R2_marginal R2_conditional 
#          0.61          19.10 


orchard_plot(mod_resp,
             mod = "ResponsePeriod",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_resp, mod = "ResponsePeriod")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
After	0.350	0.183	0.517	0.000	-1.463	2.163
During	0.494	0.317	0.671	0.000	-1.321	2.308
After-During	0.144	-0.083	0.371	0.214	-1.676	1.964

Control type

Code

# control type
# Mix only N = 1 so remove 

dat %>% filter(ControlType != "mix") ->
  data_cont

VCVs4 <- vcalc(vi = data_cont$Vd,
             cluster = data_cont$SubjectID,
             rho = 0.5)

mod_cont <- rma.mv(yi = SMD, 
               V = VCVs4, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ControlType - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = data_cont)
summary(mod_cont)

Multivariate Meta-Analysis Model (k = 639; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-891.3557  1782.7115  1794.7115  1821.4427  1794.8450   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0080  0.0893     90     no  FocalSpL 
sigma^2.2  0.1504  0.3879    115     no     RecNo 
sigma^2.3  0.6901  0.8307    639     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 636) = 13.3306, p-val < .0001

Model Results:

                        estimate      se    tval   df    pval   ci.lb   ci.ub 
ControlTypeBlank          0.4424  0.0947  4.6709  636  <.0001  0.2564  0.6283 
ControlTypedisturbance    0.4170  0.1452  2.8710  636  0.0042  0.1318  0.7022 
ControlTypeNonPred        0.3895  0.0810  4.8106  636  <.0001  0.2305  0.5485 
                            
ControlTypeBlank        *** 
ControlTypedisturbance   ** 
ControlTypeNonPred      *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_cont)*100, 2)

   R2_marginal R2_conditional 
          0.07          18.73 

Code

# R2_marginal R2_conditional 
#           0.07          18.73 

orchard_plot(mod_cont,
             mod = "ControlType",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_cont, mod = "ControlType")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Blank	0.442	0.256	0.628	0.000	-1.376	2.261
Disturbance	0.417	0.132	0.702	0.004	-1.414	2.248
NonPred	0.389	0.231	0.548	0.000	-1.426	2.205
Blank-Disturbance	-0.025	-0.346	0.295	0.876	-1.862	1.812
Blank-NonPred	-0.053	-0.260	0.155	0.617	-1.874	1.768
Disturbance-NonPred	-0.027	-0.346	0.291	0.865	-1.864	1.809

Age

Code

# age
#A = Adult, N = Nestling, J = Juveniles (including fledglings and first year birds), E = Egg

# "" only N = 2 (study = 1) so remove 

dat %>% filter(Age != "") ->
  dat_age

VCVs5 <- vcalc(vi = dat_age$Vd,
             cluster = dat_age$SubjectID,
             rho = 0.5)


# age
mod_age <- rma.mv(yi = SMD, 
               V = VCVs5, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Age - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat_age)
summary(mod_age)

Multivariate Meta-Analysis Model (k = 638; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-888.1737  1776.3474  1790.3474  1821.5117  1790.5263   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0084  0.0917     90     no  FocalSpL 
sigma^2.2  0.1498  0.3870    116     no     RecNo 
sigma^2.3  0.6914  0.8315    638     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 634) = 9.9438, p-val < .0001

Model Results:

      estimate      se    tval   df    pval    ci.lb   ci.ub      
AgeA    0.4206  0.0677  6.2136  634  <.0001   0.2877  0.5536  *** 
AgeE    0.3228  0.3452  0.9353  634  0.3500  -0.3550  1.0006      
AgeJ    0.2881  0.3484  0.8269  634  0.4086  -0.3961  0.9724      
AgeN    0.3265  0.1911  1.7087  634  0.0880  -0.0487  0.7018    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_age)*100, 2)

   R2_marginal R2_conditional 
          0.15          18.74 

Code

   # R2_marginal R2_conditional 
   #        0.15          18.74 

orchard_plot(mod_age,
             mod = "Age",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_age, mod = "Age")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.421	0.288	0.554	0.000	-1.394	2.235
E	0.323	-0.355	1.001	0.350	-1.610	2.256
J	0.288	-0.396	0.972	0.409	-1.647	2.223
N	0.327	-0.049	0.702	0.088	-1.522	2.175
A-E	-0.098	-0.782	0.587	0.779	-2.033	1.837
A-J	-0.133	-0.823	0.558	0.706	-2.070	1.805
A-N	-0.094	-0.472	0.283	0.625	-1.943	1.755
E-J	-0.035	-0.997	0.927	0.944	-2.085	2.015
E-N	0.004	-0.766	0.774	0.992	-1.963	1.971
J-N	0.038	-0.738	0.815	0.923	-1.931	2.008

Predactor type

Code

# type of predator

dat$PredTo[dat$PredTo == ""] <- NA



mod_predt <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredTo - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_predt)

Multivariate Meta-Analysis Model (k = 621; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-861.2315  1722.4629  1734.4629  1761.0219  1734.6004   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     85     no  FocalSpL 
sigma^2.2  0.1572  0.3965    110     no     RecNo 
sigma^2.3  0.6805  0.8249    621     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 618) = 12.3559, p-val < .0001

Model Results:

         estimate      se    tval   df    pval    ci.lb   ci.ub      
PredToA    0.4284  0.0814  5.2628  618  <.0001   0.2686  0.5883  *** 
PredToB    0.3772  0.3065  1.2306  618  0.2190  -0.2247  0.9791      
PredToN    0.3337  0.0915  3.6466  618  0.0003   0.1540  0.5135  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_predt)*100, 2)

   R2_marginal R2_conditional 
          0.25          18.97 

Code

# R2_marginal R2_conditional 
#           0.25          18.97 


orchard_plot(mod_predt,
             mod = "PredTo",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# result table
all_models(mod_predt, mod = "PredTo")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.428	0.269	0.588	0.000	-1.376	2.233
B	0.377	-0.225	0.979	0.219	-1.518	2.273
N	0.334	0.154	0.513	0.000	-1.473	2.140
A-B	-0.051	-0.674	0.572	0.872	-1.953	1.851
A-N	-0.095	-0.314	0.125	0.397	-1.905	1.716
B-N	-0.043	-0.672	0.585	0.892	-1.947	1.861

8 Meta-regression: multi-moderator

Code

#######################
# Mulit-variable models
#######################

dat_short$sln_duration <- scale(log(dat_short$duration_days))

mod_full <- rma.mv(yi = SMD, 
               V = VCVs, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ sln_duration*Type +
                         sln_duration*Treat_mod +
                         Sex,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat_short)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-826.1989  1652.3977  1682.3977  1748.0486  1683.2369   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     81     no  FocalSpL 
sigma^2.2  0.1540  0.3925    104     no     RecNo 
sigma^2.3  0.7169  0.8467    600     no    Obs_ID 

Test of Moderators (coefficients 2:12):
F(df1 = 11, df2 = 588) = 2.0578, p-val = 0.0216

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
intrcpt                         0.4853  0.1076   4.5122  588  <.0001   0.2741 
sln_duration                   -0.1333  0.0896  -1.4886  588  0.1371  -0.3093 
TypeLifeHistory                -0.0267  0.4683  -0.0569  588  0.9546  -0.9463 
TypePhysiology                 -0.5076  0.1840  -2.7580  588  0.0060  -0.8691 
Treat_modV                      0.0477  0.1461   0.3264  588  0.7442  -0.2393 
Treat_modAV                     0.1566  0.1656   0.9455  588  0.3448  -0.1687 
SexF                           -0.2748  0.1295  -2.1219  588  0.0343  -0.5291 
SexM                           -0.0870  0.1531  -0.5679  588  0.5703  -0.3877 
sln_duration:TypeLifeHistory    0.1502  0.4474   0.3358  588  0.7372  -0.7284 
sln_duration:TypePhysiology     0.2185  0.1995   1.0952  588  0.2739  -0.1733 
sln_duration:Treat_modV        -0.0308  0.1418  -0.2168  588  0.8284  -0.3093 
sln_duration:Treat_modAV       -0.1650  0.1804  -0.9150  588  0.3606  -0.5193 
                                ci.ub      
intrcpt                        0.6966  *** 
sln_duration                   0.0426      
TypeLifeHistory                0.8930      
TypePhysiology                -0.1461   ** 
Treat_modV                     0.3347      
Treat_modAV                    0.4818      
SexF                          -0.0205    * 
SexM                           0.2138      
sln_duration:TypeLifeHistory   1.0289      
sln_duration:TypePhysiology    0.6103      
sln_duration:Treat_modV        0.2478      
sln_duration:Treat_modAV       0.1892      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_full)*100, 2)

   R2_marginal R2_conditional 
          7.06          23.50 

Code

orchard_plot(mod_full,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

orchard_plot(mod_full,
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# interaction plot of duration with trait types
int_type <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                                   by = "Type")

bubble_plot(int_type, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days) (z-transformed)",
                     legend.pos = "top.left", condition.nrow = 3, g = TRUE)

Code

# interaction plot of duration with trait types
int_trt <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                        by = "Treat_mod")

bubble_plot(int_trt, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3, g = TRUE)

Code

# mulit-model selection
candidates <- dredge(mod_full, trace = 2)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%

Code

# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 5) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 5)) 

# relative importance of each predictor for all the models
importance <- sw(candidates)

9 Publication bias

Funnel plot: uni-moderator

Code

funnel(mod_ma2, 
       yaxis="seinv",
       # = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 16),
       xlim = c(-10,15),
       digits=c(0,1)
       )

Funnel plot: multi-moderator

Code

funnel(mod_full, 
       yaxis="seinv",
       #type = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 1.5),
       xlim = c(-5,10),
       digits=c(0,1)
       )

Egger regression: uni-moderator

Code

# Egger

dat$effectN <- (dat$Ncontrol * dat$NTreat) / (dat$Ncontrol + dat$NTreat)
dat$sqeffectN <- sqrt(dat$effectN)

mod_mae <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ sqeffectN,
                random = list(
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod_mae)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.0018  1786.0037  1796.0037  1818.2954  1796.0986   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     90     no  FocalSpL 
sigma^2.2  0.1483  0.3850    116     no     RecNo 
sigma^2.3  0.6929  0.8324    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 1.0650, p-val = 0.3025

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt      0.5044  0.1194   4.2234  638  <.0001   0.2699  0.7390  *** 
sqeffectN   -0.0291  0.0282  -1.0320  638  0.3025  -0.0844  0.0262      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod_mae,
            mod = "sqeffectN",
            group = "RecNo",
            xlab = "Effective N",
            g = TRUE)

Decline effect: uni-moderator

Code

# decline effect
mod_mad <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ Year,
                random = list(
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod_mad)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.5966  1785.1932  1795.1932  1817.4849  1795.2881   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1450  0.3808    116     no     RecNo 
sigma^2.3  0.6886  0.8298    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 1.5062, p-val = 0.2202

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt   23.4292  18.7537   1.2493  638  0.2120  -13.3973  60.2556    
Year      -0.0114   0.0093  -1.2273  638  0.2202   -0.0297   0.0069    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod_mad,
            mod = "Year",
            group = "RecNo",
            xlab = "Publication year",
            g = TRUE)

All together

Code

# full model
dat_short$effectN <- (dat_short$Ncontrol * dat_short$NTreat) / (dat_short$Ncontrol + dat_short$NTreat)
dat_short$sqeffectN <- sqrt(dat_short$effectN)

mod_fulle <- rma.mv(yi = SMD, 
                   V = VCVs, 
                   random = list(~1|FocalSpL , 
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mods =  ~ sqeffectN +
                     Year +
                     sln_duration*Type +
                     sln_duration*Treat_mod +
                     Sex,
                   test = "t",
                   method = "REML",
                   sparse = TRUE,
                   data = dat_short)
summary(mod_fulle)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-823.4875  1646.9750  1680.9750  1755.3215  1682.0525   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     81     no  FocalSpL 
sigma^2.2  0.1587  0.3983    104     no     RecNo 
sigma^2.3  0.7183  0.8475    600     no    Obs_ID 

Test of Moderators (coefficients 2:14):
F(df1 = 13, df2 = 586) = 1.7639, p-val = 0.0454

Model Results:

                              estimate       se     tval   df    pval     ci.lb 
intrcpt                        14.6850  22.7802   0.6446  586  0.5194  -30.0557 
sqeffectN                      -0.0007   0.0328  -0.0216  586  0.9827   -0.0651 
Year                           -0.0070   0.0113  -0.6221  586  0.5341   -0.0293 
sln_duration                   -0.1252   0.0936  -1.3370  586  0.1817   -0.3090 
TypeLifeHistory                -0.0173   0.4714  -0.0366  586  0.9708   -0.9431 
TypePhysiology                 -0.5322   0.1894  -2.8094  586  0.0051   -0.9043 
Treat_modV                      0.0368   0.1478   0.2489  586  0.8035   -0.2535 
Treat_modAV                     0.1375   0.1690   0.8136  586  0.4162   -0.1944 
SexF                           -0.2740   0.1300  -2.1076  586  0.0355   -0.5293 
SexM                           -0.0953   0.1543  -0.6178  586  0.5369   -0.3983 
sln_duration:TypeLifeHistory    0.1450   0.4510   0.3216  586  0.7479   -0.7407 
sln_duration:TypePhysiology     0.2438   0.2055   1.1863  586  0.2360   -0.1599 
sln_duration:Treat_modV        -0.0261   0.1440  -0.1813  586  0.8562   -0.3088 
sln_duration:Treat_modAV       -0.1497   0.1834  -0.8161  586  0.4148   -0.5099 
                                ci.ub     
intrcpt                       59.4257     
sqeffectN                      0.0637     
Year                           0.0152     
sln_duration                   0.0587     
TypeLifeHistory                0.9086     
TypePhysiology                -0.1602  ** 
Treat_modV                     0.3271     
Treat_modAV                    0.4693     
SexF                          -0.0187   * 
SexM                           0.2077     
sln_duration:TypeLifeHistory   1.0307     
sln_duration:TypePhysiology    0.6475     
sln_duration:Treat_modV        0.2566     
sln_duration:Treat_modAV       0.2105     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# small-study
bubble_plot(mod_fulle, group = "RecNo", mod = "sqeffectN", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3, g = TRUE)

Code

# decline
bubble_plot(mod_fulle, group = "RecNo", mod = "Year", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3, g = TRUE)

10 Figure gallery

Figure: alluvial plots

Figure: phylogeny + world map

Figure: main and different modality

Figure: different moderators

11 R Session Informtion

Code

# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.3.1 (2023-06-16)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: pacman(v.0.5.1), pander(v.0.6.5), kableExtra(v.1.3.4), MuMIn(v.1.47.5), emmeans(v.1.8.7), clubSandwich(v.0.5.10), ape(v.5.7-1), easyalluvial(v.0.3.1), ggalluvial(v.0.12.5), alluvial(v.0.1-2), patchwork(v.1.1.3), metafor(v.4.2-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), orchaRd(v.2.0), lme4(v.1.1-34), Matrix(v.1.6-1), here(v.1.0.1), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.2), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.3) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): gridExtra(v.2.3), sandwich(v.3.0-2), rlang(v.1.1.1), magrittr(v.2.0.3), ggridges(v.0.5.4), compiler(v.4.3.1), mgcv(v.1.8-42), systemfonts(v.1.0.4), vctrs(v.0.6.3), rvest(v.1.0.3), pkgconfig(v.2.0.3), fastmap(v.1.1.1), labeling(v.0.4.2), utf8(v.1.2.3), rmarkdown(v.2.24), prodlim(v.2023.03.31), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), nloptr(v.2.0.3), xfun(v.0.40), jsonlite(v.1.8.7), recipes(v.1.0.7), highr(v.0.10), parallel(v.4.3.1), R6(v.2.5.1), stringi(v.1.7.12), RColorBrewer(v.1.1-3), parallelly(v.1.36.0), boot(v.1.3-28.1), rpart(v.4.1.19), estimability(v.1.4.1), Rcpp(v.1.0.11), knitr(v.1.43), future.apply(v.1.11.0), zoo(v.1.8-12), splines(v.4.3.1), nnet(v.7.3-19), timechange(v.0.2.0), tidyselect(v.1.2.0), rstudioapi(v.0.15.0), yaml(v.2.3.7), timeDate(v.4022.108), codetools(v.0.2-19), listenv(v.0.9.0), lattice(v.0.21-8), withr(v.2.5.0), coda(v.0.19-4), evaluate(v.0.21), future(v.1.33.0), survival(v.3.5-5), xml2(v.1.3.5), pillar(v.1.9.0), stats4(v.4.3.1), generics(v.0.1.3), rprojroot(v.2.0.3), mathjaxr(v.1.6-0), hms(v.1.1.3), munsell(v.0.5.0), scales(v.1.2.1), minqa(v.1.2.5), globals(v.0.16.2), xtable(v.1.8-4), class(v.7.3-22), glue(v.1.6.2), tools(v.4.3.1), data.table(v.1.14.8), webshot(v.0.5.5), gower(v.1.0.1), mvtnorm(v.1.2-2), grid(v.4.3.1), ipred(v.0.9-14), colorspace(v.2.1-0), nlme(v.3.1-162), beeswarm(v.0.4.0), vipor(v.0.4.5), latex2exp(v.0.9.6), cli(v.3.6.1), fansi(v.1.0.4), viridisLite(v.0.4.2), svglite(v.2.1.1), lava(v.1.7.2.1), gtable(v.0.3.3), digest(v.0.6.33), progressr(v.0.14.0), farver(v.2.1.1), htmlwidgets(v.1.6.2), htmltools(v.0.5.6), lifecycle(v.1.0.3), hardhat(v.1.3.0), httr(v.1.4.7) and MASS(v.7.3-60)